
import sys
import time
import readgenome

KNOWN_GENE_TABLE = "/home/bowtie/Bowtie_Illumina_Alignment/references/knownGene.table.txt"

resComp = { "A":"T", "T":"A", "C":"G", "G":"C", "N":"N", "a":"t", "t":"a", "g":"c", "c":"g", "n":"n" }

def revComp(st):
  strev = ""
  for i in xrange(0, len(st)):
    strev += resComp[st[len(st)-i-1]]
  return strev

if __name__=="__main__":
  t0 = time.time()
  numGenes = 0
  allChrs = {}
  for chr in ["chr%s"%x for x in (range(1,23) + ["X", "Y"])]:
    allChrs[chr] = True
  for line in open(KNOWN_GENE_TABLE):
    if line.startswith("#"): continue
    cols = line.strip().split("\t")
    (name, chrom, strand, txStart, txEnd, cdsStart, cdsEnd, exonCount, exonStarts, exonEnds, proteinID, alignID) = tuple(cols)
    if not allChrs.has_key(chrom):
      print >> sys.stderr, "# Skipping %s on unknown reference %s" % (name, chrom)
      continue
    txStart = long(txStart)
    txEnd = long(txEnd)
    exonStarts = [long(x) for x in exonStarts.split(",")[:-1]]
    exonEnds = [long(x) for x in exonEnds.split(",")[:-1]]
    exons = [(exonStarts[i], exonEnds[i]) for i in xrange(int(exonCount))]
    fullSeq = readgenome.readGenomicSequence((chrom, txStart+1, txEnd+1))
    transcript = ""
    rng = ""
    for (spos, epos) in exons:
      transcript += fullSeq[spos-txStart:epos-txStart]
      if rng=="":
        rng += chrom + ":%d-%d" % (spos+1, epos)
      else:
        rng += ",%d-%d" % (spos+1, epos)
    if strand=="-":
      transcript = revComp(transcript)
    print ">" + name + " range=" + rng + " strand=" + strand
    pos = 0
    while pos<len(transcript):
      print transcript[pos:pos+50]
      pos+=50
    numGenes += 1
    if 0==numGenes%100:
      print >> sys.stderr, "# Processed %d genes in time %f" % (numGenes, time.time() - t0)
